Skip to content

Conversation

@sanguinariojoe
Copy link
Collaborator

Related:

#290
#261

The new flexible states are here! With these new states it is rather easy to install new variables, which can have fixed lengths (real, vec, vec6...) or variable lengths depending on the instance (e.g. lines as vec3 and rods as a combination of vec7 and vec6).

It is fairly optimized, avoiding some copies here and there. However, the implementation comes to a performance cost, becoming as twice slower as the previous implementation. I suppose some other optimizations can be made to avoid some copies. I just do not know if it is worth the effort, as long as we are bringing --hopefully soon-- a time scheme which will make simulations waaaaaay faster.

@sanguinariojoe sanguinariojoe added this to the 2.4.0 milestone Feb 24, 2025
@sanguinariojoe sanguinariojoe self-assigned this Feb 24, 2025
@sanguinariojoe
Copy link
Collaborator Author

Windows compiling is failing... I will fix that tomorow

@sanguinariojoe sanguinariojoe force-pushed the flexstate branch 3 times, most recently from a0e96f3 to fa23498 Compare February 25, 2025 08:36
source/Line.cpp Outdated
}

void
Line::setState(const StateVarRef pos, const StateVarRef vel)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this meant to replace Line::setState above? If so then we will need to make similar changes to Bodies, Rods, and Points I assume

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is done:

https://github.com/core-marine-dev/MoorDyn/blob/flexstate/source/Body.hpp#L383
https://github.com/core-marine-dev/MoorDyn/blob/flexstate/source/Rod.hpp#L473

On those cases it is more explicit though. It is totally possible to use list to make it more flexible anyway.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, so in the header files it converts these new state variables to the old inputs to the functions? If feels redundant to have two Line::setState functions

source/Time.cpp Outdated
Update(0.0, 0);
CalcStateDeriv(0);
r[0] = r[0] + rd[0] * dt;
r(0)->get<list>("pos") += dt * rd(0)->get<list>("vel");
Copy link
Collaborator

@RyanDavies19 RyanDavies19 Feb 25, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just so I have it clear, if we were to add an additional state for lines, would it live in r(i)->get<list>("pos") where pos is now a list with more elements or would it be its own line here like r(0)->get<list>(<new state name>) += dt * rd(0)->get<list>(<new state deriv name);

If it is the latter then we would need to do it for every integration scheme it seems

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Both are possible. I like more to have a new state var name, and if so create a set of macros to integrate things... Something like:

#define STATE_VAR_DECL(i)  \
    STATE_VAR_BLOCK(T) r ## i = r(i)->get<list>("pos"); \
    STATE_VAR_BLOCK(T) u ## i = r(i)->get<list>("vel"); \
    STATE_VAR_BLOCK(T) l ## i = r(i)->get<list>("length");

#define STATE_DERIV_DECL(T, i)  \
    STATE_VAR_BLOCK(T) drdt ## i = rd(i)->get<list>("vel"); \
    STATE_VAR_BLOCK(T) dudt ## i = rd(i)->get<list>("acc"); \
    STATE_VAR_BLOCK(T) dldt ## i = rd(i)->get<list>("dldt");

#define STATE_INTEGRATE_INPLACE(ivar, ideriv, dt)  \
    r ## ivar += dt * drdt ## ideriv;  \
    u ## ivar += dt * dudt ## ideriv;  \
    l ## ivar += dt * dldt ## ideriv;

#define STATE_INTEGRATE(ivar_out, ivar_in, ideriv, dt)  \
    r ## ivar_out = r ## ivar_in + dt * drdt ## ideriv;  \
    u ## ivar_out = u ## ivar_in + dt * dudt ## ideriv;  \
    l ## ivar_out = l ## ivar_in + dt * dldt ## ideriv;

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Said that, packing everything on a single state variable --here it will set the list length as 6 to hold both pos and vel-- would make things way faster (no more ::get<>() and less copies).

Slightly less idiomatic maybe, but again I think it can be fixed with some macros.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually the latter way is pretty sexy as well. Let me know the way you wanna go please (and do not worry if I have to rewrite all, I am happily doing that)

Copy link
Collaborator

@RyanDavies19 RyanDavies19 Feb 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the latter might be a better bet, provided I am understanding it right. When setting up the line state, you would do something like

virtual void AddLine(Line* obj, bool viv, bool visco)
	{
		try {
			Scheme::AddLine(obj);
		} catch (...) {
			throw;
		}
		// Build up the states and state derivatives
		AS_STATE(_r[0])->addLine(obj);
		AS_STATE(_r[0])->setListLength("pos", 3, obj);
		AS_STATE(_r[0])->setListLength("vel", 3, obj);
                if (viv) AS_STATE(_r[0])->setListLength("phase", 1, obj);
                if (visco) AS_STATE(_r[0])->setListLength("length", 1, obj);

And then in the time schemes we could just integrate r? Still not super clear on what that would look like but happy to try.

What I am trying to get at is: how do we do this so we only allocate and integrate these extra line states if the user requests that feature? In the MD-F code I do this:

Set state vector size: https://github.com/RyanDavies19/MoorDynF_ryan/blob/1c883b541da4f91c84cd4d072dd8cc8926852ad9/modules/moordyn/src/MoorDyn.f90#L1489
Allocate line variables if needed: https://github.com/RyanDavies19/MoorDynF_ryan/blob/1c883b541da4f91c84cd4d072dd8cc8926852ad9/modules/moordyn/src/MoorDyn_Line.f90#L168

@RyanDavies19
Copy link
Collaborator

@sanguinariojoe, on first glance this looks pretty nice, thanks again for putting it together. I left two clarifying comments above that would be good to discuss. I will take a stab at getting the VIV and viscoelastic work functioning with this approach and report back (likely with some more questions).

@RyanDavies19
Copy link
Collaborator

Also, pinging @AlexWKinley if you have bandwidth to review that would be great because this changes a lot of what you did on #269

@AlexWKinley
Copy link
Contributor

I appreciate the ping.
I've looked through this a little bit.
Code wise, there's nothing here I'm particularly opposed to. What I am more opposed to is the significant performance penalty.

I haven't run any benchmarks myself, or done any profiling, so I'm speaking with no authority on what the root of the performance regression is.

One thing that stands out is the use of std::map. std::map is generally implemented as a red black tree, which gives good theoretical performance, but at the cost of practical performance (std::map is slow in the same way a linked list is slow).

Given that the set of possible state variables is known, I think you should be able to get away with just a std::vector, and maybe an enum for readability.

A rough outline of how that might look is

// In state.hpp replace:

/// The map of available variables
std::map<std::string, VarBase*> vars;

/// The map of available variable types
std::map<std::string, VarBase::types> types;

// with:

/// List of variable names
std::vector<std::string> varNames;

/// List of variables, order matches varNames
std::vector<VarBase*> vars;

/// List of variable types
std::vector<VarBase::types> types;

And then you'd make the required changes elsewhere, such that something like this from Time.cpp could be changed

// You could define these enums somewhere
enum StateVars {
    POS = 0,
    VEL = 1,

};
enum DStateVars {
    VEL = 0,
    ACC = 1,
}

// then in RK2::Step for example, you could have
auto r0 = r(0)->get<list>(StateVars::POS);
auto u0 = r(0)->get<list>(StateVars::VEL);
auto r1 = r(1)->get<list>(StateVars::POS);
auto u1 = r(1)->get<list>(StateVars::VEL);
auto drdt0 = rd(0)->get<list>(DStateVars::VEL);
auto dudt0 = rd(0)->get<list>(DStateVars::ACC);
auto drdt1 = rd(1)->get<list>(DStateVars::VEL);
auto dudt1 = rd(1)->get<list>(DStateVars::ACC);

That all being said, I'll try and get some time today to run some benchmarks and do some profiling, to see if there are any easy improvements to be made.

@sanguinariojoe
Copy link
Collaborator Author

sanguinariojoe commented Feb 26, 2025 via email

@sanguinariojoe
Copy link
Collaborator Author

I just want to let you know that I have ready the new system. I shall still test and fix it.

It is taking a bit longer than expected because one of my kids is sick. Hopefully tomorrow I am uploading it.

@RyanDavies19
Copy link
Collaborator

I just want to let you know that I have ready the new system. I shall still test and fix it.

It is taking a bit longer than expected because one of my kids is sick. Hopefully tomorrow I am uploading it.

Hopefully they feel better soon!

@sanguinariojoe
Copy link
Collaborator Author

sanguinariojoe commented Mar 5, 2025 via email

@sanguinariojoe
Copy link
Collaborator Author

Here we go!

@RyanDavies19 , to create new state variables you can use the methods:

moordyn::Line::stateN()
moordyn::Line::stateDims()

For instance, to add a the length of each segments you can edit those methods as:

/** @brief Get the number of state variables required by this instance
	 * @return The number of internal nodes
	 * @warning This function shall be called after ::setup()
	 */
	inline const size_t stateN() const { return getN(); }

	/** @brief Get the dimension of the state variable
	 * @return 3 components for positions and 3 components for velocities, i.e.
	 * 6 components
	 * @note In future developments this function might return different
	 * numbers depending on the line configuration. e.g. the viscoelasticity.
	 * @warning This function shall be called after ::setup()
	 */
	inline const size_t stateDims() const { return 7; }

Then I suppose you must handle that on

https://github.com/core-marine-dev/MoorDyn/blob/flexstate/source/Line.cpp#L682
https://github.com/core-marine-dev/MoorDyn/blob/flexstate/source/Line.cpp#L780

@RyanDavies19
Copy link
Collaborator

@sanguinariojoe, at a very quick glance this looks really nice and easy to work with! Love that there is a single state and state derivative object passed into things we can access from. Keeps all the nice features here while also using a similar structure as MD-F. I will take a closer look and play around with this more. I'll let you know here if I have any questions.

@RyanDavies19
Copy link
Collaborator

RyanDavies19 commented Mar 10, 2025

Okay @sanguinariojoe, I was working with this and I think we are quite close. For implementing the VIV model, the existing structure works and I can put it together pretty easily. Working on that now.

For the viscoelastic model, we need a state for every segment, rather than every internal node. This means that Nstates needs to be different for this dimension. For example, in this approach a normal line has a N-1 x 6 instance (where N is the number of nodes). Disregarding the VIV model, if we want to include the viscoelastic model, then we need to have an additional N states, i.e. a N x 1 array. This raises an issue because either we need to make a N x 7 instance and have a couple unused spaces in the instance (and strange indexing for the segments), or a slightly revised approach. What are your thoughts on adding a Nx1 instance when using the viscoelastic model, so that we have an instance for the Line Nodes and an instance for the Line Segments? We would have to figure out how this ties to the line object of course, but it would allow for future developments that want to expand the states on the segments (for example, if someone implemented a generalized Maxwell model for viscoelastic behavior). This would probably mean modifying the AddLine and RemoveLine functions in Time.hpp to also handle the Segment states. If users weren't using the viscoelastic model, we could easily disable any unnecessary instances. Let me know if this approach sounds good to you.

A side question, I noticed that states are being allocated and assigned to all objects in the system including coupled objects. Do we need to do this? MD-F does not add coupled object states to the state list. I know the performance improvements will probably be minimal, but it would probably help. I would think it would be as simple as changing the calls to t_integrator.Add<object> to iterate over the free objects. I.e. change:

	t_integrator.SetGround(GroundBody);
	for (auto obj : BodyList)
		t_integrator.AddBody(obj);
	for (auto obj : RodList)
		t_integrator.AddRod(obj);
	for (auto obj : PointList)
		t_integrator.AddPoint(obj);
	for (auto obj : LineList)
		t_integrator.AddLine(obj);

to the following (in icLegacy and icStationary in MoorDyn2.cpp):

	t_integrator.SetGround(GroundBody);
	for (unsigned int i : FreeBodyIs)
		t_integrator.AddBody(BodyList[i]);
	for (unsigned int i : FreeRodIs)
		t_integrator.AddRod(RodList[i]);
	for (unsigned int i : FreePointIs)
		t_integrator.AddPoint(PointList[i]);
	for (auto obj : LineList)
		t_integrator.AddLine(obj);

@sanguinariojoe
Copy link
Collaborator Author

For the viscoelastic model, we need a state for every segment, rather than every internal node. This means that Nstates needs to be different for this dimension. For example, in this approach a normal line has a N-1 x 6 instance (where N is the number of nodes). Disregarding the VIV model, if we want to include the viscoelastic model, then we need to have an additional N states, i.e. a N x 1 array. This raises an issue because either we need to make a N x 7 instance and have a couple unused spaces in the instance (and strange indexing for the segments), or a slightly revised approach. What are your thoughts on adding a Nx1 instance when using the viscoelastic model, so that we have an instance for the Line Nodes and an instance for the Line Segments? We would have to figure out how this ties to the line object of course, but it would allow for future developments that want to expand the states on the segments (for example, if someone implemented a generalized Maxwell model for viscoelastic behavior). This would probably mean modifying the AddLine and RemoveLine functions in Time.hpp to also handle the Segment states. If users weren't using the viscoelastic model, we could easily disable any unnecessary instances. Let me know if this approach sounds good to you.

In my mind a N x 7 is the way to go. It is very true that we have some unused slots on the line state matrix, but it is just a chunk on the last row of each line. Regarding the indexing, it works actually the same way it does now inside Line.cpp. I see no difference.

The main merit of just growing stateN() by one --when viscoelasticity is enabled-- is the simplicity then, i.e. the states are still matrices, and they work the same fashion for lines, points, rods and bodies. No need for special cases.
On top of that, they still enjoy the benefits of the stateDims() function, so lines with viscoelastic models --eventually including a generalized Maxwell one-- may ask for a N x M state matrix, becoming both N and M different for each line.

On the other hand, if those new states does not have a derivative (I do not know how they are computed, so you tell me), it would be totally worthy to summon another state variable, so we'll have sort of kinematics state variable and constants state variable:

https://github.com/core-marine-dev/MoorDyn/blob/flexstate/source/Time.hpp#L792

	/// The list of kinematic states
	std::array<moordyn::state::State*, NSTATE> _r;

	/// The list of kinematic state derivatives
	std::array<moordyn::state::State*, NDERIV> _rd;

	/// The list of constants
	std::array<moordyn::state::State*, NDERIV> _c;

Let me know please

A side question, I noticed that states are being allocated and assigned to all objects in the system including coupled objects. Do we need to do this? MD-F does not add coupled object states to the state list. I know the performance improvements will probably be minimal, but it would probably help. I would think it would be as simple as changing the calls to t_integrator.Add<object> to iterate over the free objects. I.e. change:

	t_integrator.SetGround(GroundBody);
	for (auto obj : BodyList)
		t_integrator.AddBody(obj);
	for (auto obj : RodList)
		t_integrator.AddRod(obj);
	for (auto obj : PointList)
		t_integrator.AddPoint(obj);
	for (auto obj : LineList)
		t_integrator.AddLine(obj);

to the following (in icLegacy and icStationary in MoorDyn2.cpp):

	t_integrator.SetGround(GroundBody);
	for (unsigned int i : FreeBodyIs)
		t_integrator.AddBody(BodyList[i]);
	for (unsigned int i : FreeRodIs)
		t_integrator.AddRod(RodList[i]);
	for (unsigned int i : FreePointIs)
		t_integrator.AddPoint(PointList[i]);
	for (auto obj : LineList)
		t_integrator.AddLine(obj);

Yeah, totally... I think we went the current way for simplicity again, but I think we can totally move to just free entities. I just wonder about interacting with moorpy... Anyway, I will give it a shot after fixing the other issue

@RyanDavies19
Copy link
Collaborator

In my mind a N x 7 is the way to go. It is very true that we have some unused slots on the line state matrix, but it is just a chunk on the last row of each line. Regarding the indexing, it works actually the same way it does now inside Line.cpp. I see no difference.

The main merit of just growing stateN() by one --when viscoelasticity is enabled-- is the simplicity then, i.e. the states are still matrices, and they work the same fashion for lines, points, rods and bodies. No need for special cases. On top of that, they still enjoy the benefits of the stateDims() function, so lines with viscoelastic models --eventually including a generalized Maxwell one-- may ask for a N x M state matrix, becoming both N and M different for each line.

On the other hand, if those new states does not have a derivative (I do not know how they are computed, so you tell me), it would be totally worthy to summon another state variable, so we'll have sort of kinematics state variable and constants state variable:

https://github.com/core-marine-dev/MoorDyn/blob/flexstate/source/Time.hpp#L792

	/// The list of kinematic states
	std::array<moordyn::state::State*, NSTATE> _r;

	/// The list of kinematic state derivatives
	std::array<moordyn::state::State*, NDERIV> _rd;

	/// The list of constants
	std::array<moordyn::state::State*, NDERIV> _c;

Let me know please

Let's go with the Nx7 approach for now then. The viscoelastic state does have a derivative, no constants needed. I do like the simplicity of it. I will work on getting both of these implemented and will update my PR. Is the best approach to merge this one and then run tests again on mine? Or is it better to open one to your flexstate branch?

Yeah, totally... I think we went the current way for simplicity again, but I think we can totally move to just free entities. I just wonder about interacting with moorpy... Anyway, I will give it a shot after fixing the other issue

Sounds good, that can be a project for further down the line. I have gotten some interest recently about exporting from MoorPy to MoorDyn for a project happening at NREL so thanks for setting it up!

@sanguinariojoe
Copy link
Collaborator Author

Yeah, totally... I think we went the current way for simplicity again, but I think we can totally move to just free entities. I just wonder about interacting with moorpy... Anyway, I will give it a shot after fixing the other issue

Sounds good, that can be a project for further down the line. I have gotten some interest recently about exporting from MoorPy to MoorDyn for a project happening at NREL so thanks for setting it up!

I gave it a fast shot, and the results were not good. A lot of tests started failing. I suppose we can address this later. Would you open an issue?

Is the best approach to merge this one and then run tests again on mine? Or is it better to open one to your flexstate branch?

This works fine, so I think it will be easier and faster to get it merged

@RyanDavies19
Copy link
Collaborator

@sanguinariojoe,

I'll open an issue about getting that done in future work.

Also, happy to report the new states work great. I've got the VIV and viscoelastic model running well and their performance has been verified. I'm going to clean some stuff up and then update my PR. This one is good to merge from my perspective!

@RyanDavies19
Copy link
Collaborator

@sanguinariojoe

One other thought, has this been made compatible with line failures? I know those have not been upkept much but was curious how this ties in there.

@RyanDavies19
Copy link
Collaborator

RyanDavies19 commented Mar 12, 2025

@sanguinariojoe, one more things that could be nice to add here, but could be left for later. It would be good to keep setState, getStateDeriv, and initialize consistent between the objects. For the changes in #295, I had lines accept a InstanceStateVarView as inputs to these functions (see https://github.com/RyanDavies19/MoorDynC_ryan/blob/d7f309d5bfc509524051d4efc35218df336eb821/source/Line.cpp#L759). I tried getting the same going with Points, Bodies, and Rods, but it gets tricky because there are multiple calls to initialize outside of the time scheme files (ex. bodies initializing attached objects, or points and rods initializing in the main Init function). I think the easiest solution would be to find a way to have those objects provide a InstanceStateVarView, using a modified version of getState. I took a stab at that, but wasn't having much luck. Before banging my head on the wall too much, I thought I would ask you if you had an idea of a simple solution.

@sanguinariojoe
Copy link
Collaborator Author

cdac02f @RyanDavies19 ?

If you confirm that's ok, I am merging tomorrow

@RyanDavies19
Copy link
Collaborator

RyanDavies19 commented Mar 17, 2025

@sanguinariojoe,

Those look good. I had been thinking it would be nice to pass the state directly into the <object>::Initialize function (in .cpp), but I understand why you set it like this. It does help make the system more simple. I think for initialization we should go with your approach.

For the <object>::SetState calls I guess the best approach is to document things well. In my mind, SetState and getStateDeriv should have identical input/output structure, because they are core to understanding the structure of the code. I have made the change for bodies and rods, points are tricky because of failures. I will save them for later when I am updating the line failures capability.

@RyanDavies19
Copy link
Collaborator

This is ready to merge once we figure out the failing build tests. Could this be related to #296? Prior to your last commit they were passing (see be0a344).

@sanguinariojoe
Copy link
Collaborator Author

This is ready to merge once we figure out the failing build tests. Could this be related to #296? Prior to your last commit they were passing (see be0a344).

Nah, it is because Rust. I just upgraded the FindRust.cmake module and everything is working alright again.

I am thus merging this!

@sanguinariojoe sanguinariojoe merged commit 63f8784 into FloatingArrayDesign:dev Mar 18, 2025
7 of 8 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants